On the Numerical Calculation of the Roots of Special Functions Satisfying Second Order Ordinary Differential Equations

We describe a method for calculating the roots of special functions satisfying second order ordinary differential equations. It exploits the recent observation that the solutions of equations of this type can be represented via nonoscillatory phase functions, even in the high-frequency regime. Our approach requires O(1) operations per root and achieves near machine precision accuracy. Moreover, despite its great generality, our approach is competitive with (and in many cases, faster than) specialized, state-of-the-art methods for the construction of Gaussian quadrature rules of large orders when it used in such a capacity. The performance of the scheme is illustrated with several numerical experiments.